4.1 Вероятностный подход в ML

Как описать привычные модели на языке статистики. Оптимизация функции потерь vs оценка максимального правоподобия

В этом разделе мы посмотрим на те же самые модели машинного обучения, но с другой стороны: будем интерпретировать их как вероятностные.

В первом параграфе мы расскажем, как обращаться с вероятностными моделями, и покажем, что привычный вам подбор параметров модели с помощью минимизации функции потерь соответствует подбору параметров методом максимального правдоподобия. Это даст возможность транслировать в мир ML известные результаты о свойствах оценок максимального правдоподобия, но в то же время и обнажит их недостатки. Благодаря этому мы сможем по-новому взглянуть на логистическую регрессию и с новым пониманием сформулировать её обобщение — generalized linear model (GLM).

По ходу дела мы обнаружим, что большинство классификаторов, хоть и делают вид, что предсказывают корректные вероятности, на самом деле вводят в заблуждение.

В третьем параграфе мы поговорим о том, как проверить отклонение предсказанных значений от истинных вероятностей и как поправить ситуацию.

Далее мы обсудим генеративный подход к классификации и разберём несколько примеров генеративных моделей, после чего перейдём к байесовскому подходу оценивания параметров, который, хоть зачастую и трудно осуществим вычислительно, однако обладает большей теоретической стройностью. Он позволяет оценивать распределение параметров и предсказаний (например, уверенность в нашей оценке), а кроме того — даёт нам возможность измерить качество модели, не прибегая к проверке на тестовой выборке.

Если вы готовы — давайте приступим!

Случайность как источник несовершенства модели

Практически любая наша модель — несовершенна. Но объяснять это несовершенство можно по-разному.

Представим, что мы решаем задачу регрессии yx,wy\simeq \langle x, w\rangle: например, пытаемся по университетским оценкам выпускника предсказать его годовую зарплату. Ясно, что точная зависимость у нас не получится как минимум потому, что мы многого не знаем о выпускнике: куда он пошёл работать, насколько он усерден, как у него с soft skills и так далее. Как же нам быть?

Первый вариант — просто признать, что мы не получим идеальную модель, но постараться выучить оптимальную, насколько это возможно. То есть приблизить таргет предсказаниями наилучшим образом с точки зрения какой-то меры близости, которую мы подберём из экспертных соображений.

Так мы получаем простой инженерный подход к машинному обучению: есть формула, в которой присутствуют некоторые параметры (ww), есть формализация того, что такое «приблизить» (функция потерь) — и мы бодро решаем задачу оптимизации по параметрам.

Второй вариант — свалить вину за неточности наших предсказаний на случайность. В самом деле: если мы что-то не можем измерить, то для нас это всё равно что случайный фактор. В постановке задачи мы заменяем приближённое равенство yx,wy\simeq \langle x, w\rangle на точное

y=(x,w, искажённое шумом ε)y=(\langle x, w\rangle, \text { искажённое шумом } \varepsilon)

Например, это может быть аддитивный шум (чаще всего так и делают):

y=x,w+εy = \langle x, w\rangle + \varepsilon

где ε\varepsilon — некоторая случайная величина, которая представляет этот самый случайный шум. Тогда получается, что для каждого конкретного объекта xix_i соответствующий ему истинный таргет — это сумма xi,w\langle x_i, w\rangle и конкретной реализации шума ε\varepsilon.

При построении такой модели мы можем выбирать различные распределения шума, кодируя тем самым, какой может быть ошибка. Чаще всего выбирают гауссовский шум: εN(0,σ2)\varepsilon\sim\mathcal{N}(0,\sigma^2) с некоторой фиксированной дисперсией σ2\sigma^2 — но могут быть и другие варианты.

Проиллюстрируем, как ведут себя данные, подчиняющиеся закону y=ax+b+εy = ax + b + \varepsilon, εN(0,σ2)\varepsilon\sim\mathcal{N}(0, \sigma^2):

9

Вопрос на подумать. Зачем человеку может прийти в голову предположить, что в модели линейной регрессии yXw+εy\sim Xw + \varepsilon шум ε\varepsilon имеет распределение Лапласа? А распределение Коши? Чем свойства таких моделей будут отличаться от свойств модели с нормальным шумом?

Ответ (не открывайте сразу; сначала подумайте сами!)

Давайте посмотрим, как выглядят плотности этих трёх распределений:

9

Распределение Лапласа имеет «более тяжёлые хвосты», чем нормальное: это значит, что плотность медленнее падает с удалением от среднего. Таким образом, этому распределению могут подчиняться данные, в которых имеются выбросы. Если не гнаться за строгостью, можно сказать, что модель с нормальным шумом будет пытаться объяснить выбросы, меняя под них ww, тогда как лапласовский шум потерпит их, не подгоняя ww.

У распределения Коши хвосты «ещё более тяжёлые», что, в теории, даёт возможность модели с таким шумом описывать даже ещё более шумные данные.

Проиллюстрируем датасеты, сгенерированные из моделей с каждым из типов шума: нормальным, лапласовским и Коши.

9

Как вы могли заметить, в каждом из подходов после того, как мы зафиксировали признаки (то есть координаты xix_i), остаётся своя степень свободы: в инженерном это выбор функции потерь, а в вероятностном — выбор распределения шума.

Дальше в этом параграфе мы увидим, что на самом деле эти два подхода глубинным образом связаны между собой, причём выбор функции потерь — это в некотором смысле то же самое, что выбор распределения шума.

Условное распределение на таргет, непрерывный случай

Допустим, что мы исследуем вероятностную модель таргета с аддитивным шумом

y=fw(x)+ε,y = f_w(x) + \varepsilon,

где fwf_w — некоторая функция, не обязательно линейная с (неизвестными пока) параметрами ww, а ε\varepsilon — случайный шум с плотностью распределения εpε(t)\varepsilon\sim p_{\varepsilon}(t). Для каждого конкретного объекта xix_i значение fw(xi)f_w(x_i) — это просто константа, но для yiy_i оно превращается в случайную величину, зависящую от xix_i (и ещё от ww, на самом деле).

Таким образом, можно говорить об условном распределении

py(yx,w)p_y(y \vert x, w)

Для каждого конкретного xix_i и ww распределение соответствующего yiy_i — это просто pε(yfw(xi))p_{\varepsilon}(y - f_{w}(x_i)), ведь yfw(X)=εy - f_w(X) = \varepsilon.

Пример. Рассмотрим вероятностную модель y=x,w+εy = \langle x, w\rangle + \varepsilon, где εN(0,σ2)\varepsilon\sim\mathcal{N}(0, \sigma^2). Тогда для фиксированного xix_i имеем yi=xi,w+εy_i = \langle x_i, w\rangle + \varepsilon. Поскольку xi,w\langle x_i, w\rangle — константа, мы получаем

yiN(xi,w,σ2).y_i\sim\mathcal{N}(\langle x_i, w\rangle, \sigma^2).

Это можно записать и так:

p(yixi,w)N(yixi,w,σ2),p(y_i\vert x_i, w)\sim\mathcal{N}(y_i\vert\langle x_i, w\rangle, \sigma^2),

где выражение справа — это значение функции плотности нормального распределения с параметрами xi,w,σ2\langle x_i, w\rangle, \sigma^2 в точке yiy_i. В частности, xi,w=E(yixi)\langle x_i, w\rangle = \mathbb{E}(y_i\vert x_i).

Более сложные вероятностные модели

На самом деле, мы можем для нашей задачи придумывать любую вероятностную модель py(yx,w)p_y(y \vert x, w), не обязательно вида y=fw(X)+εy = f_w(X) + \varepsilon.

Представьте, что мы хотим предсказывать точку в плоскости штанг, в которую попадает мячом бьющий по воротам футболист. Можно предположить, что она имеет нормальное распределение со средним (цель удара), которое определяется ситуацией на поле и состянием игрока, и некоторой дисперсией (то есть скалярной ковариационной матрицей), которая тоже зависит от состояния игрока и ещё разных сложных факторов, которые мы объявим случайными.

Состояние игрока — это сложное понятие, но, вероятно, мы можем выразить его, зная пульс, давление и другие физические показатели. В свою очередь, ситуацию на поле можно описать, как функцию от позиций и движений других игроков, судьи и зрителей — но всего не перечислишь, поэтому нам снова придётся привлекать случайность. Таким образом, мы получаем то, что называется графической моделью:

9

Здесь стрелки означают статистические зависимости, а отсутствие стрелок — допущение о статистической независимости. Конечно же, это лишь допущение, принятое нами для ограничения сложности модели: ведь пульс человека и давление взаимосвязаны, равно как и поведение различных игроков на поле. Но мы уже обсуждали, что каждая модель, в том числе и вероятностная, является лишь приблизительным отражением бесконечно сложного мира. Впрочем, если у нас много вычислительных ресурсов, то никто не мешает нам попробовать учесть и все пропущенные сейчас зависимости.

Расписав всё по определению условной вероятности, мы получаем следующую вероятностную модель:

9

в которой, конечно же, мы должны все вероятности расписать через какие-то понятные и логически обоснованные распределения — но пока воздержимся от этого.

Оценка максимального правдоподобия = оптимизация функции потерь

Мы хотим подобрать такие значения параметров ww, для которых модель py(yx,w)p_y(y \vert x, w) была бы наиболее адекватна обучающим данным. Суть метода максимального правдоподобия (maximum likelihood estimation) состоит в том, чтобы найти такое ww, для которого вероятность (а в данном, непрерывном, случае плотность вероятности) появления выборки y={y1,,yN}y = \{y_1, \ldots, y_N\} была бы максимальной, то есть

w^MLE=argmaxwp(yX,w)\widehat{w}_{MLE} = \underset{w}{\operatorname{argmax}}p(y \vert X, w)

Величина p(yX,w)p(y \vert X, w) называется функцией правдоподобия (likelihood). Если мы считаем, что все объекты независимы, то функция правдоподобия распадается в произведение:

p(yX,w)=p(y1x1,w)p(yixi,w)p(y \vert X, w) = p(y_1 \vert x_1, w) \cdot\ldots\cdot p(y_i \vert x_i, w)

Теперь, поскольку перемножать сложно, а складывать легко (и ещё поскольку мы надеемся, что раз наши объекты всё-таки наблюдаются в природе, их правдоподобие отлично от нуля), мы переходим к логарифму функции правдоподобия:

l(yX,w)=logp(y1x1,w)++logp(yixi,w)l(y \vert X,w) = \log{p(y_1 \vert x_1, w)} + \ldots + \log{p(y_i \vert x_i, w)}

эту функцию мы так или иначе максимизируем по ww, находя оценку максимального правдоподобия w^\hat{w}.

Как мы уже обсуждали выше, p(yixi,w)=pε(yfw(xi))p(y_i \vert x_i, w) = p_{\varepsilon}(y - f_{w}(x_i)), то есть

l(yX,w)=i=1Nlogpε(yifw(xi))l(y \vert X,w) = \sum\limits_{i=1}^N\log{p_{\varepsilon}(y_i - f_w(x_i))}

Максимизация функции правдоподобия соответствует минимизации

i=1N[logpε(yifw(xi))]\sum\limits_{i=1}^N\left[-\log{p_{\varepsilon}(y_i - f_w(x_i))}\right]

а это выражение можно интерпретировать, как функцию потерь. Вот и оказывается, что подбор параметров вероятностей модели с помощью метода максимального правдоподобия — это то же самое, что «инженерная» оптимизация функции потерь. Давайте посмотрим, как это выглядит в нескольких простых случаях.

Пример. Давайте предположим, что наш таргет связан с данными вот так:

yi=xi,w+εy_i = \langle x_i, w \rangle + \varepsilon

где εN(0,σ2)\varepsilon\sim\mathcal{N}(0, \sigma^2), то есть

p(ε)=12πσ2exp(ε22σ2)p(\varepsilon) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{\varepsilon^2}{2\sigma^2}\right)

Случайная величина yiy_i получается из шума ε\varepsilon сдвигом на постоянный вектор xi,w\langle x_i, w \rangle, так что она тоже распределена нормально с той же дисперсией σ2\sigma^2 и со средним xi,w\langle x_i, w \rangle

p(yixi,w)=12πσ2exp((yixi,w)22σ2)p(y_i\vert \langle x_i, w \rangle) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{(y_i - \langle x_i, w \rangle)^2}{2\sigma^2}\right)

Правдоподобие выборки имеет вид

p(yX,w)=i=1Np(yixi,w)=i=1N12πσ2exp((yiw,xi)22σ2)p(y\vert X, w) = \prod_{i=1}^N p(y_i \vert x_i, w) = \prod_{i=1}^N \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left(-\frac{(y_i-\langle w,x_i\rangle)^2}{2\sigma^2}\right)

Логарифм правдоподобия можно переписать в виде

l(yX,w)=i=1N(log(2πσ2)(yiw,xi)22σ2)l(y \vert X,w) = \sum_{i=1}^N \left(-\log({\sqrt{2 \pi \sigma^2}}) -\frac{(y_i-\langle w,x_i\rangle)^2}{2\sigma^2}\right)

Постоянными слагаемыми можно пренебречь, и тогда оказывается, что максимизация этой величины равносильна минимизации

i=1N(yiw,xi)2 \sum_{i=1}^N (y_i-\langle w,x_i\rangle)^2

Мы получили обычную квадратичную функцию потерь. Итак, обучать вероятностную модель линейной регрессии с нормальным шумом — это то же самое, что учить «инженерную» модель с функцией потерь MSE.

Вопрос на подумать. Какая вероятностная модель соответствует обучению линейной регрессии с функцией потерь MAE

i=1Nyiw,xi?\sum_{i=1}^N \vert y_i-\langle w,x_i\rangle\vert?

Ответ (не открывайте сразу; сначала подумайте сами!)

Минимизация функции потерь MAE соответствует максимизации

i=1N[yiw,xi]\sum_{i=1}^N\left[-\vert y_i-\langle w,x_i\rangle\vert\right]

Мы хотим найти такое распределение, для которого эта штука является с точностью до константы логарифмом функции правдоподобия. Что ж, возьмём экспоненту:

exp[i=1Nyiw,xi]=i=1Nexp(yiw,xi)\text{exp}\left[-\sum_{i=1}^N\vert y_i-\langle w,x_i\rangle\vert\right] = \prod_{i=1}^N\text{exp}\left(-\vert y_i-\langle w,x_i\rangle\vert\right)

Если теперь это умножить на (12)N\left(\frac12\right)^{N}, то мы получим функцию правдоподобия для распределения Лапласа:

i=1N12exp(yiw,xi)=i=1NLaplace(yiw,xi)\prod_{i=1}^N\frac12\text{exp}\left(-\vert y_i-\langle w,x_i\rangle\vert\right) = \prod_{i=1}^NLaplace\left(y_i-\langle w,x_i\rangle\right)

Итак, учить «инженерную» модель с функцией потерь MAE — это то же самое, что обучать вероятностную модель линейной регрессии с лапласовским шумом.

Предсказание в вероятностных моделях

Теперь представим, что параметры подобраны, и подумаем о том, как же теперь делать предсказания.

Рассмотрим модель линейной регрессии

y=x,w+ε,εN(0,σ2)y = \langle x, w\rangle + \varepsilon,\quad\varepsilon\sim\mathcal{N}(0,\sigma^2)

Если ww известен, то для нового объекта x0x_0 соответствующий таргет имеет вид

y0=x0,w+εN(x0,w,σ2)y_0 = \langle x_0, w\rangle + \varepsilon\sim\mathcal{N}(\langle x_0, w\rangle, \sigma^2)

Таким образом, y0y_0 дан нам не точно, а в виде распределения (и логично: ведь мы оговорились выше, что ответы у нас искажены погрешностью, проинтерпретированной, как нормальный шум). Но что делать, если требуют назвать конкретное число? Кажется логичным выдать условное матожидание E(y0x0)=x0,w\mathbb{E}(y_0\vert x_0) = \langle x_0, w\rangle, тем более что оно совпадает с условной медианой и условной модой этого распределения.

Если же медиана, мода и математическое ожидание различаются, то можно выбрать что-то из них с учётом особенностей задачи. Но на практике в схеме yf(x)+εy\sim f(x) + \varepsilon чаще всего рассматривают именно симметричные распределения с нулевым матожиданием, потому что для них f(x)f(x) совпадает с условным матожиданием E(yx)\mathbb{E}(y\vert x) и является логичным точечным предсказанием.

Приведём пример. Допустим шум ε\varepsilon был бы из экспоненциального распределения. Тогда f(x)f(x) была бы условным минимумом распределения. В принципе, можно придумать задачу, для которой такая постановка (предсказание минимума) была бы логичной. Но это всё же довольно экзотическая ситуация. Приводим для сравнения модели с нормальным, лапласовским и экспоненциальным шумом:

9

Условное распределение на таргет, дискретный случай

Допустим, мы имеем дело с задачей классификации с KK классами. Как мы можем её решать? Самый наивный вариант — научиться по каждому объекту xix_i предсказывать некоторое число для каждого класса, и у кого число больше — тот класс и выбираем! Наверное, так можно сделать, если мы придумаем хорошую функцию потерь. Но сразу в голову приходит мысль: почему бы не начать предсказывать не просто число, а вероятность?

Таким образом, задача классификации сводится к предсказанию

P(yi=kxi)P(y_i = k \vert x_i)

и как будто бы выбору класса с наибольшей вероятностью. Впрочем, как мы увидим дальше, всё не всегда работает так просто.

Одну такую модель — правда, только для бинарной классификации — вы уже знаете. Это логистическая регрессия:

P(yi=1xi,w)=11+exi,w,P(yi=0xi,w)=e(xi,w)1+exi,w=11+exi,wP(y_i = 1 \vert x_i,w) = \frac{1}{1+e^{-\langle x_i, w\rangle}},\quad P(y_i = 0 \vert x_i,w) = \frac{e^{-(x_i, w)}}{1+e^{-\langle x_i, w\rangle}} = \frac{1}{1+e^{\langle x_i, w\rangle}}

которую также можно записать в виде

yixiBern(11+exi,w)y_i \vert x_i \sim \color{red}{Bern}\left(\frac{1}{1+e^{-\langle x_i, w\rangle}}\right)

где Bern(p)\color{red}{Bern}(p) — распределение Бернулли с параметром pp.

Нахождение вероятностей классов можно разделить на два этапа:

 Находим xi логиты (xi,w,xi,w)σ(σ(xi,w),σ(xi,w))\begin{aligned} &\text { Находим }\\ &x_i \xrightarrow{\text { логиты }}\left(-\left\langle x_i, w\right\rangle,\left\langle x_i, w\right\rangle\right) \xrightarrow{\sigma}\left(\sigma\left(-\left\langle x_i, w\right\rangle\right), \sigma\left(\left\langle x_i, w\right\rangle\right)\right) \end{aligned}

где, напомним, σ\sigma — это сигмоида:

σ(t)=11+et\sigma(t) = \frac{1}{1+e^{-t}}

Сигмоида тут не просто так. Она обладает теми счастливыми свойствами, что

  • монотонно возрастает;

  • отображает всю числовую прямую на интервал (0,1)(0,1);

  • σ(x)=1σ(x)\sigma(-x) = 1 - \sigma(x).

Вот такой вид имеет её график:

9

Иными словами, с помощью сигмоиды можно делать «вероятности» из чего угодно, то есть более или менее для любого отображения fwf_w (из признакового пространства в R\mathbb{R}) с параметрами ww построить модель бинарной классификации:

P(yi=0xi,w)=σ(fw(xi)),P(yi=1xi,w)=σ(fw(xi)).P(y_i = 0 \vert x_i, w) = \sigma(f_w(-x_i)),\quad P(y_i = 1 \vert x_i, w) = \sigma(f_w(x_i)).

Как и в случае логистической регрессии, такая модель равносильна утверждению о том, что

fw(xi)=logp(y=1xi,w)p(y=0xi,w).f_w(x_i) = \log{\frac{p(y = 1 \vert x_i,w)}{p(y = 0 \vert x_i, w)}}.

Похожим способом можно строить и модели для многоклассовой классификации. В этом нам поможет обобщение сигмоиды, которое называется softmax:

softmax(t1,,tK)=(et1k=1Ketk,,etKk=1Ketk)softmax(t_1,\ldots,t_K) = \left(\frac{e^{t_1}}{\sum_{k=1}^Ke^{t_k}},\ldots,\frac{e^{t_K}}{\sum_{k=1}^Ke^{t_k}}\right)

А именно, для любого отображения fwf_w из пространства признаков в RK\mathbb{R}^K мы можем взять модель

(P(yi=kxi,w))k=1K=softmax(fw(xi))\left(P(y_i = k \vert x_i, w)\right)^K_{k=1} = softmax(f_w(x_i))

Если все наши признаки — вещественные числа, а fw(xi)=xiWf_w(x_i) = x_iW — просто линейное отображение, то мы получаем однослойную нейронную сеть

(P(yi=kxi,w))k=1K=softmax(xiW)\left(P(y_i = k \vert x_i, w)\right)^K_{k=1} = softmax(x_iW)

9

Предостережение. Всё то, что мы описали выше, вполне работает на практике (собственно, классификационные нейросети зачастую так и устроены), но корректным не является.

В самом деле, мы говорим, что строим оценки вероятностей P(yi=kxi,w)P(y_i = k \vert x_i, w), но для подбора параметров используем не эмпирические вероятности, а только лишь значения argmaxk P(yi=kxi,w)\underset{k}{\operatorname{argmax}} \ P(y_i = k \vert x_i, w), то есть метки предсказываемых классов. Таким образом, при обучении мы не будем различать следующие две ситуации:

9

Это говорит нам о некоторой неполноценности такого подхода.

Заметим ещё вот что. В случае бинарной классификации выбор предсказываемого класса как argmaxkP(yi=kxi,w)\underset{k}{\operatorname{argmax}} P(y_i=k \vert x_i,w) равносилен выбору того класса, для которого P(yi=kxi,w)>12P(y_i=k \vert x_i,w) > \frac{1}{2}. Но если наши оценки вероятностей неадекватны, то этот вариант проваливается, и мы встаём перед проблемой выбора порога: каким должно быть значение t^\widehat{t}, чтобы мы могли приписать класс 1 тем объектам xix_i, для которых σ(fw(xi))>t^\sigma(f_w(x_i)) > \widehat{t}?

В одном из следующих параграфов мы обсудим, как всё-таки правильно предсказывать вероятности.

Чтобы добавить в заметки выделенный текст, нажмите Command + E
Предыдущий параграф3.4. Подбор гиперпараметров

Как эффективно подбирать значения гиперпараметров модели и не переобучиться при этом

Следующий параграф4.2. Экспоненциальный класс распределений и принцип максимальной энтропии

Самые главные семейства распределений в жизни любого data scientist’а

Вступайте в сообщество хендбука

Здесь можно найти единомышленников, экспертов и просто интересных собеседников. А ещё — получить помощь или поделиться знаниями.